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1. Introduction 

Image deconvolution has been an active research field for several decades and recent contributions 
can be found in papers such as |[I]-[3l. Examples of application are medical imaging, astronomy, 
nondestructive testing and more generally imagery problems. In these applications, degradations 
induced by the observation instrument limit the data resolution while the need of precise interpre- 
tation can be of major importance. For example, this is particularly critical for long-wavelength 
astronomy (see e.g., JH). In addition, the development of a high quality instrumentation system 
must rationally be completed by an equivalent level of quality in the development of data process- 
ing methods. Moreover, even for poor performance systems, the restoration method can be used to 
bypass instrument limitations. 
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When the deconvolution problem is ill-posed a possible solution relies on regularization, i.e., in- 
troduction of information in addition to the data and the acquisition model jSMJ. As a consequence 
of regularization, deconvolution methods are specific to the class of image in accordance with the 
introduced information. From this standpoint, the present paper is dedicated to relatively smooth 
images encountered for numerous applications in imagery [|4l|7l[8l. The second order consequence 
of ill-posedness and regularization is the need to balance the compromise between different sources 
of information. 

In the Bayesian approach EIHl, information about unknowns is introduced by means of proba- 
bilistic models. Once these models are designed, the next step is to build the a posteriori law, given 
the measured data. The solution is then defined as a representative point of this law and the two 
most classical are (1) the maximizer, and (2) the mean. From a computational standpoint, the first 
leads to a numerical optimization problem and the latter leads to a numerical integration problem. 
However, the resulting estimate depends on two sets of variables in addition to the data. 

1. Firstly, the estimate naturally depends on the response of the instrument at work, namely the 
point spread function (PSF). The literature is predominantly devoted to deconvolution in the 
case of known PSF. On the contrary, the present paper is devoted to the case of unknown 
or poorly known PSF and there are two main strategies to tackle its estimation from the 
available data set (without extra measurements). 

(i) In most practical cases, the instrument can be modeled using physical operating descrip- 
tion. It is thus possible to find the equation for the PSF, at least in a first approximation. 
This equation is usually driven by a relatively small number of parameters. It is a com- 
mon case in optical imaging where a Gaussian-shaped PSF is often used ifTOl . It is also 
the case in other fields: interferometry [il 111 , magnetic resonance force microscopy [[T2ll . 
fluorescence microscopy lfT3l .. . . Nevertheless, in real experiments, the parameter values 
are unknown or imperfectly known and need to be estimated or adjusted in addition to 
the image of interest: the question is namely myopic deconvolution. 

(ii) The second strategy forbears the use of the parametric PSF deduced from the physical 
analysis and the PSF then naturally appears in a non-parametric form. Practically, the 
non-parametric PSF is unknown or imperfectly known and needs to be estimated in 
addition to the image of interest: the question is referred to as blind deconvolution for 
example in interferometry in"4] - [T7l . 

From an inference point of view, the difficulty of both myopic and blind problems lies in 
the possible lack of information resulting in ambiguity between image and PSF, even in the 
noiseless case. In order to resolve the ambiguity, information must be added |l3l[l8l and it 
is crucial to make inquiries based on any available source of information. To this end, the 
knowledge of the parametric PSF represents a precious means to structure the problem and 
possibly resolve the degeneracies. Moreover, due to instrument design process, a nominal 
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value as well as an uncertainty are usually available for the PSF parameters. 
In addition, from a practical and algorithmic standpoint, the myopic case, i.e., the case of 
parametric PSF, is often more difficult due to the non-linear dependence of the observation 
model with respect to the PSF parameters. On the contrary, the blind case, i.e., the case of 
non-parametric PSF, yields a simpler practical and algorithmic problem since the observation 
model remains linear w.r.t. the unknown elements given the object. 

Despite the superior technical difficulty, the present paper is devoted to the myopic format 
since it is expected to be more efficient than the blind format from an information standpoint. 
Moreover, the blind case has been extensively studied and a large amount of paper is available 
lfT9l - l2T| . while the myopic case has been less investigated, though it is of major importance. 
2. Secondly, the solution depends on the probability law parameters named hyperparameters 
(means, variances, parameters of correlation matrix,. . . ). These parameters adjust the shape 
of the laws and in the same time they tune the compromise between the information provided 
by the a priori and the information provided by the data. In real experiments, their values are 
unknown and need to be estimated: the question is namely unsupervised deconvolution. 
For both families of parameters (PSF parameters and hyperparameters), two approaches are 
available. In the first one, the parameter values are empirically tuned or estimated in a preliminary 
step (with Maximum Likelihood Ol or calibration [|22| for example), then the values are used in a 
second step devoted to image restoration given the parameters. In the second one, the parameters 
and the object are jointly estimated [|2l[T9l. 

For the myopic problem, Jalobeanu et al. [|23l address the case of a symmetric Gaussian PSF. 
The width parameter and the noise variance are estimated in a preliminary step by Maximum- 
Likelihood. A recent paper [[24| addresses the estimation of a Gaussian blur parameter, as in our 
experiment, with an empirical method. They found the Gaussian blur parameter by minimizing the 
absolute derivatives of the restored images Laplacian. 

The present paper addresses the myopic and unsupervised deconvolution problem. We propose 
a new method that jointly estimates the PSF parameters, the hyperparameters, and the image of 
interest. It is built in a coherent and global framework based on an extended a posteriori law 
for all the unknown variables. The posterior law is obtained via the Bayes rule, founded on a 
priori laws: Gaussian for image and noise, uniform for PSF parameters and gamma or Jeffreys for 
hyperparameters . 

Regarding the image prior law, we have paid special attention to the parametrization of the 
covariance matrix in order to facilitate law manipulations such as integration, conditioning or hy- 
perparameter estimation. The possible degeneracy of the a posteriori law in some limit cases is 
also studied. 

The estimate is chosen as the mean of the posterior law and is computed using Monte-Carlo 
simulations. To this end, Monte-Carlo Markov chain (MCMC) algorithms (25] enable to draw 
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samples from the posterior distribution despite its complexity and especially the non-linear depen- 
dence w.r.t. the PSF parameters. 

The paper is structured in the following manner. Sec. [21 presents the notations and states the 
problem. The three following sections describe our methodology: firstly the Bayesian probabilistic 
models are detailed in Sec. [31 then a proper posterior law is established in Sec. [H an MCMC 
algorithm to compute the estimate is described in Sec. [51 Numerical results are shown in Sec. [6l 
Finally, Sec.[71is devoted to conclusion and perspectives. 

2. Notations and convolution model 

Consider N pixels real square images represented in lexicographic order by vector x E R^, with 
generic elements x„. The forward model writes 

y = H^x + e (1) 

where y E is the vector of data, a convolution matrix, x the image of interest and e the 
modelization errors or the noise. Vector w E R^ stands for the PSF parameters, such as width or 
orientation of a Gaussian PSF. 

The matrix is block-circulant with circulant-block (BCCB) for computational efficiency 
of the convolution in the Fourier space. The diagonalization |[26ll of writes Ah = FH^F'^ 
where F is the unitary Fourier matrix and f is the transpose conjugate symbol. The convolution, 
in the Fourier space, is then 

y = KHX + °e (2) 

where x = Fx, y = Fy and e = Fe are the 2D discrete Fourier transform (dft-2d) of image, 
data and noise, respectively. 

Since Ah is diagonal, the convolution is computed with a term-wise product in the Fourier 
space. There is a strict equivalence between a description in spatial domain (Eq. Qi) and in Fourier 
domain (Eq. ©). Consequently, for coherent description and computational efficiency, all the de- 
velopments are equally done in the spatial space or in the Fourier space. 

For notational convenience, let us introduce the component at null-frequency G R and the 
vector of component at non-null frequencies x^ E C^~^ so that the whole set of components writes 

o ro o T 

X = [Xo,X^.\. 

Let us note 1 the vector of N components equal to l/N, so that t^x is the empirical mean level 

o o o 

of the image. The Fourier components are the 1„ and we have: Iq = 1 and 1„ = for n 0. 
Moreover, Ai = FU^F"^ is a diagonal matrix with only one non-null coefficient at null frequency. 

3. Bayesian probabilistic model 

This section presents the prior law for each set of parameters. Regarding the image of interest, in 
order to account for smoothness, the law introduces high-frequency penalization through a differ- 
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ential operator on the pixel. A conjugate law is proposed for the hyperparameters and a uniform 
law is considered for the PSF parameters. 

Moreover, we have paid a special attention to the image prior law parametrization. In the next 
section we present several parametrization in order to facilitate law manipulations such as inte- 
gration, conditioning or hyperparameter estimation. Moreover, the correlation matrix of the image 
law may become singular in some limit cases resulting in a degenerated prior law (when p{x) = 
for all X G R.^). Based on this parametrization. Sec. |4] studies the degeneracy of the posterior in 
relation with the parameters of the prior law. 

3.A. Image prior law 

The probability law for the image is a Gaussian field with a given precision matrix P parametrized 
by a vector 7. The pdf reads 

1 



x^Px 
2 



(3) 



For computational efficiency, the precision matrix is designed (or approximated) in a toroidal man- 
ner, and it is diagonal in the Fourier domain Ap = FPF^ . Thus, the law for x also writes 



p{x\^) = (27r)-^/2^et[F]det[Ap]i/Met[F^] exp 
= (27r)-^/Met[Ap]i/2 exp 



--x^F^ApFx 
2 



1 o t . o 

— X Apx 
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(4) 
(5) 



and it is sometimes referred to [27 J as a Whittle approximation (see also 11281 p. 133]) for the 
Gaussian law. The filter obtained for fixed hyperparameters is also the Wiener- Hunt filter ll29l . as 
described in Sec. 15.Al 

This paper focuses on smooth images, thus on positive correlation between pixels. It is intro- 
duced by high-frequencies penalty using any circulant differential operator: p-th differences be- 
tween pixels, Laplacian, Sobel. . . The differential operator is denoted by D and its diagonalized 
form by Ao = FDF^ . Then, the precision matrix writes P = ^iD^D and its Fourier counterpart 
writes 

Ap = -^fiA^jjAo = diag ( 0,7i|rfip, . . • ,7iMAf-in (6) 



where 71 is a positive scale factor, diag builds a diagonal matrix from elementary components and 

o 

dn is the n-th DFT-2D coefficient of D. 

Under this parametrization of P, the first eigenvalue is equal to zero corresponding to the ab- 
sence of penalty for the null frequency xq, i.e., no information accounted for about the empirical 
mean level of the image. As a consequence, the determinant vanishes det[P] = resulting in a 
degenerated prior. To manage this difficulty, several approaches have been proposed. 

Some authors ElllOl still use this prior despite its degeneracy and this approach can be analyzed 
in two ways. 
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1 . On the one hand, it can be seen as a non-degenerated law for , the set of non-null frequency 
components only. In this format, the prior does not affect any probability to the null frequency 
component and the Bayes rule does not apply to this component. Thus, this strategy yields 
an incomplete posterior law, since the null frequency is not embedded in the methodology. 

2. On the other hand, it can be seen as a degenerated prior for the whole set of frequencies. The 
application of the Bayes rule is then somewhat confusing due to degeneracy. In this format, 
the posterior law cannot be guaranteed to remain non-degenerated. 

Anyway, none of the two standpoints yields a posterior law that is both non-degenerated and 
addressing the whole set of frequencies. 

An alternative parametrization relies on the energy of x. An extra term 79/, tuned by 70 > 0, in 
the precision matrix OTI . introduces information for all the frequencies including xq. The precision 
matrix writes 

Ap = lol + 71 AI^Ad 

= diag (70,70 + 7i Ml P,- • • , 7o + 7iM7v-in (7) 



with a determinant 

N-l , s 

det[Ap] = n 7o + 7iMnr • (8) 

n=0 ^ ' 

The obtained Gaussian prior is not degenerated and undoubtedly leads to a proper posterior. Nev- 
ertheless, the determinant Eq. ([8]) is not separable in 70 and 71. Consequently, the conditional 
posterior for these parameters is not a classical law and future development will be more difficult. 
Moreover, the non-null frequencies controlled by two parameters 70 and 71 

p(a;|7o,7i) = p(a;o|7o)p(a?*|7o,7i)- (9) 

The proposed approach to manage the degeneracy relies on the addition of a term for the null 
frequency only Aj = diag (1, 0, . . . , 0) 

Ap = 7oA1Ai + 71AI5AC. (10) 
= diag (7o,7i|c?i|^ . . . ,7i|rfAr_i|^ 



The determinant has a separable expression 



N-\ 

det[Ap]=7o7f-' n (11) 



n=l 



i.e., the precision parameters have been factorized. In addition, each parameter controls a different 
set of frequencies: 

p(a;|7o,7i) =p(a;o|7o)p(a;*|7i) , 
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7o drives the empirical mean level of the image Xq and 71 drives the smoothness x^, of the image. 
With the Fourier precision structure of Eq. (flOl) . we have the non-degenerated prior law for the 
image that addresses separately all the frequencies with a factorized partition function w.r.t. (70, 71) 



p{x\7o,li) = (2vr) 



-N/2 



n=l 



70||o ,|2 

-ylkoll 



Tin . o 

2 II -'^-X?*''^* 



(12) 



where Au^ is obtained from Ad without the first line and column. The next step is to write the a 
priori law for the noise in an explicit form and the other parameters, including the law parameters 
7 and the instrument parameters w. 



3.B. Noise and data laws 

From a methodological standpoint, any statistic can be included for errors (measurement and 
model errors). It is possible to account for correlations in the error process or to account for a 
non-Gaussian law, e.g., Laplacian law, generalized Gaussian law, or other laws based on robust 
norm,. . . In the present paper, the noise is modeled as zero-mean white Gaussian vector with un- 
known precision parameter 7^ 



(27r)-^/Sf/'exp 



Consequently, the likelihood for the parameters given the observed data writes 



p{y\x,%,w) = (27r) 



— \\v 

2 



(13) 



(14) 



It naturally depends on the image x, on the noise parameter 7^ and on the PSF parameters w 
embedded in Hw. It clearly involves a least squares discrepancy that can be rewritten in the Fourier 
domain: \\y — H^xW^ = \\y — Ajfa;||^. 



3.C. Hyperparameters law 

A classical choice for hyperparameter law relies on conjugate prior [|32l : the conditional posterior 
for the hyperparameters is in the same family as its prior. It results in practical and algorithmic 
facilities: update of the laws amounts to update of a small number of parameters. 

The three parameters 70, 71 and 7^ are precision parameters of Gaussian laws Eq. (fT2l) and (fT4l) 
and a conjugate law for these parameters is the Gamma law (see Appendix |B]). Given parameters 
(a,, /3j), for z = 0, 1 or e, the pdf reads 



ai — l 

:li exp 



-7,/A),V7. G [0,+oo[ 



(15) 



In addition to computational efficiency, the law allows for non-informative priors. With specific 
parameter values, one obtains two improper non-informative prior : the Jeffreys' law ^(7) = 1/7 
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and the uniform law ^(7) = W[o,+oo[(7) with (a,, set to (0, +00) and (1, +00), respectively. 
Jeffreys' law is a classical law for the precisions and is considered as non-informative ll33l . This 
law is also invariant to power transformations: the law of 7" ll33l[34l is also a Jeffreys' law. For 
these reasons development is done using the Jeffreys' law. 

3.D. PSF parameters law 

Regarding the PSF parameters w, we consider that the instrument design process or a physical 
study provides a nominal value w with uncertainty 6, that is to say w G [itJ — (5 , w + 8\. The 
"Principle of Insufficient Reason" [|33l leads to a uniform prior on this interval 

'p{w) =U^^s{w) (16) 

where Uw,6 is a uniform pdf on [Uf — 5 , w + 8\. Nevertheless, within the proposed framework, 
the choice is not limited and other laws, such as Gaussian, are possible. Anyway other choices do 
not allow easier computation because of the non-linear dependency of the observation model w.r.t. 
PSF parameters. 



4. Proper posterior law 



At this point, the prior law of each parameter is available: the PSF parameters, the hyperparameters 
and the image. Thus, the joint law for all the parameters is built by multiplying the likelihood 
Eq. ([HI) and the a priori laws Eq. ([HI), and ([HI) 



p{x, 7„ 70, 71, w, y) = p{y\x, 7,, w)p{x\-fo, li)p{le)p{l(>)v{li)v{'w) 



and explicitly 

p(a;,7„7o,7i,'u;,^) 



(27r)-^nl-/Mn 



/3f^r(«,)/3o""r(«o)/3rr(«i: 



% _ To _ 7i 



exp 



7e||?, A o ||2 70||o ||2 7l 



— \\y-AHx\\ - —\\xo\\ - yW-^dx 



^l|2 



According to the Bayes rule, the a posteriori law reads 

p{x,'ye,7o,7u'w,y) 



where p{y) is a normalization constant 



piy) 



p{y) = J p{y,x, J, w) dxdj dw. 



(17) 



(18) 



(19) 



(20) 
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As described before, setting 70 = leads to degenerated prior and joint laws. However, when 
the observation system preserves the null frequency 70 can be considered as a nuisance parameter. 
In addition, only prior information on the smoothness is available. 

In Bayesian framework, a solution to eliminate the nuisance parameters is to integrate them out 
in the a posteriori law. According to our parametrization Sec. 13. Al the integration of 70 is the 
integration of a Gamma law. Application of Appendix IB. B I on 70 in the a posteriori law Eq. (fT9l ) 
provides 

.0 ,0. p{xo)p{y,x^,%,-fi,w\xo) 
P{x, 7o 7i, w\y) = — — — (21) 

/ p{xo)p{y, a?*, 7e, 7i, wlxo) d% d7i dw dx^ dxo 

with 



p{xo) = I pixo\lo)p{lo) d-fo 

o o2\ -"0-1/2 

1 + ^ . (22) 



Now the parameter is integrated, the parameters and (3o are set to remove the null frequency 
penalization. Since we have > and /3o > we get (1 + f3Qxl/2)~'^°~^^^ < 1 and the joint law 
is majored 

1 + -y^ piy,x^,%,'yi,w\°Xo) <piy,x^,'y^,'yi,w\°Xo). (23) 

Consequently, by the dominated convergence theorem Il35ll . the limit of the law with ao — )■ 1 
and (3q can be placed under the integral sign at the denominator. Then the null-frequency 
penalization p{xo) from the numerator and denominator are removed. It is equivalent with the 
integration of the 70 parameter under a Dirac (see appendix |B]). The equation is simplified and the 
integration with respect to a;o in the denominator Eq. (|20|) 



/ p{y\x,%,w)p{x^\-fi)p{-fi,-f^,w)dxo (X / p{yo\xo,-fe,w)dxo (24) 

2 



oc / exp 

/R 



dxo (25) 



converges if and only if /iq 7^ 0: the null frequency is observed. If this condition is met, Eq. (121)) 
with /3o = and ao = 1 is a proper posterior law for the image, the precision parameters and the 
PSF parameters. In other words, if the average is observed, the degeneracy of the a priori law is 
not transmitted to the a posteriori law. 
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Then, the obtained a posteriori law writes 
p[x,-i,,-ii,w\y) = — 

p{y) 

oc 7r^^/'~Sr^"'^''"'^/'"'W^,.(«^) (26) 

Finally, inference is done on this law Eq. (|26|) . If the null frequency is not observed, or information 
must be added, the previous Eq. ( fT9l ) can be used. 





7e 


7i 


exp 







5. Posterior mean estimator and law exploration 

This section presents the algorithm to explore the posterior law Eq. ( fT9l ) or (l26l) and to compute an 
estimate of the parameters. For this purpose, Monte Carlo Markov chain is used to provide samples. 
Firstly, the obtained samples are used to compute different moments of the law. Afterwards, they 
are also used to approximate marginal laws as histograms. These two representations are helpful 
to analyse the a posteriori law, the structure of the available information and the uncertainty. They 
are used in Sec. l6.C.2l to illustrate the mark of the ambiguity in the myopic problem. 

Here, the samples of the a posteriori law are obtained by a Gibbs sampler ll25l[36l[37l : it con- 
sists in iteratively sampling the conditional posterior law for a set of parameters given the other 
parameters (obtained at previous iteration). Typically, the sampled laws are the law of x, 7j and 
w. After a burn-in time, the complete set of samples are under the joint a posteriori law. The three 
next sections present each sampling step. 

5.A. Sampling the image 

The conditional posterior law of the image is a Gaussian law 

&('=-^^)~p(^|^,7f,7f,7f\^('^0 (27) 
~Ar(/^('=+i),E('=+^)). (28) 

The covariance matrix is diagonal and writes 

S(^+i) = (7f lAg)^ + 7riAir + 7fVz.lT' (29) 



and the mean 



^(fc+i)=^W5.(/c+i)^t^W^_ (30) 



where f is the transpose conjugate symbol. The vector fi^'^'^^^ is the regularized least square solu- 
tion at the current iteration (or the Wiener-Hunt filter). Clearly, if the null-frequency is not observed 

o 

ho = and if 70 = 0, the covariance matrix E is not invertible and the estimate is not defined as 
described Sec. IH 
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Finally, since the matrix is diagonal, the sample x is obtained by a term-wise product of 
Fe (where e is white Gaussian) with the standard deviation matrix ('S^'^'^^A ^ followed by the 



addition of the mean /x*^^"^^^ also computed with term-wise products Eq. (l30l) . Consequently, the 
sampling of the image is effective even with high-dimensional object. 

5.B. Sampling precision parameters 

The conditional posterior laws of the precisions are Gamma corresponding to their prior law with 
parameters updated by the likelihood 

ir'^-vU\yrx^"^"\^^'^] (31) 



-Q{iMr'\Pt'')- (32) 

For 7e, 7o and 71 the parameters law are, respectively, 

ap+i) = a, + N/2 and /3f = (/^^^ + \\\y - Ag^a^^'+'V) , (33) 

a^'' = ao + 1/2 and = (/^o"^ + (5:^^')')" , (34) 

a['^'^ = a, + {N- l)/2 and /^f = (/Jf^ + ^||Az.5.^'+'V)"' • (35) 
In the case of Jeffreys' prior, the parameters are 

^ik+i) ^ and = 2/\\y - Ag^a^^'+'V, (36) 

«r'^ = l/2 and P^'^'^ = 2/{xi''-'^)\ (37) 

= (iV - l)/2 and = 2/||Aoaj^'+'V. (38) 



Remark 1 — If the a posteriori law Eq. A26v without 70 is considered, there is no need to sample 
this parameter (Eq. and 07l) are not useful) and 7o'''' = in Eq. (|29l). 

5. C. Sample PSF parameters 

The conditional law for PSF parameters writes 

«;(^+i)~pfic||/,x^'+'\7f+'^) (39) 

(40) 



oc exp 



2 

where parameters w are embedded in the PSF Ah- This law is not standard and intricate: no algo- 
rithm exists for direct sampling and we use the Metropolis-Hastings (M.-H.) method to bypass this 
difficulty. In M.-H. algorithm, a sample Wp is proposed and accepted with a certain probability. 
This probability depends on the ratio between the likelihood of the proposed value and the likeli- 
hood of the current value w^'^\ In practice, in the independent form described in appendix O with 
prior law as proposition law, it is divided in several steps. 
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1. Proposition: Sample a proposition 

Wp p{w) =U[ab]{w). (41) 

2. Probability of acceptation: Calculate the criterion 

J {^^'\ ^p) = ^ (ll^ - Ah,.(^) ^^'^'V - 11^ - Ah,., x^'^'^r) . (42) 

3. Update: Sample t ~ U[o ij and takes 

I > otherwise. 

5.D. Empirical mean 

The sampling of x, 7 and it? are repeated iteratively until the law has been sufficiently explored. 



These samples 



follow the global a posteriori law of Eq. (1191 ). By the large num 



bers law, the estimate, defined as the posterior mean, is approximated by 



X 



1 K-l , ■ 

K to . 



(44) 



As described by Eq. (I44|) . to obtain an estimate of the image in the spatial space, all the computation 
are achieved recursively in the Fourier space with a single IFFT at the end. An implementation 
example in pseudo code is described Fig.|9l 

6. Deconvolution results 

This section presents numerical results obtained by the proposed method. In order to completely 
evaluate the method, true value of all parameters x, w, % but also 71,70 is needed. In order to 
achieve this, an entirely simulated case is studied: image and noise are simulated under their re- 
spective prior laws Eq. (fT2l) and (fT3l) with given values of 70, 71 and 7^. Thanks to this protocol, 
all experimental conditions are controlled and the estimation method is entirely evaluated. 

The method has also been applied in different conditions (lower signal to noise ratio, broader 
PSF, different and realistic (non- simulated) images, . . . ) and showed similar behaviour. However, 
in the case of realistic images, since the true value of the hyperparameters 70 and 71 is unknown, 
the evaluation cannot be complete. 

6.A. Practical experimental conditions 

Concretely, a 128 x 128 image is generated in the Fourier space as the product of a complex 
white Gaussian noise and the a priori standard deviation matrix S = (7oAi + 7iA^A£,)~^^^, 
given by Eq. (fTOl) . The chosen matrix Xd results from the fft-2d of the Laplacian operator 
[0 1 0; 1 — 4 1; 1 0] /8 and the parameter values are 70 = 1 and 71 = 2. 
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These parameters provide the image shown in Fig. |l(a)| : it is an image with smooth features 
similar to a cloud. Pixels have numerical values between —100 and 150, and the profile line 68 
shows fluctuations around a value of —40. 

The a priori law for the hyperparameters are set to the non-informative Jeffreys' law by fixing 
the (ai, (3i) to (0, +oo), as explained in Sec. I3.CI In addition, the PSF is obtained in the Fourier 
space by discretization of a normalized Gaussian shape 

^(t'a, 1^13) = exp ^ — 27r^ ^^^{wa cos^ (f + sin^ (p) 

+ Vpiwa sin^ ip + wp cos^ ip) 

+ 2i/q,i//3 sin cosy? (wq, — w^) (45) 

with frequencies (z/q, 1^/3) G [—0.5; 0.5]^. This low-pass filter, illustrated in Fig.[2l is controlled by 
three parameters: 

• two width parameters Wa and Wjs set to 20 and 7, respectively. Their a priori laws are uniform: 
p{wa) = W[i9 2i]{wa) and p{wi3) = U[q 8](^a) Corresponding to an uncertainty of about 5% 
and 15% around the nominal value (see Sec lB.DI) . 

• a rotation parameter ip set to 7r/3. The a priori law is also uniform p{ip) = U[jr/4 n/2]i^) 
corresponding to 50% uncertainty. 

Then, the convolution is computed in the Fourier space and the data are obtained by adding 
a white Gaussian noise with precision 7^ = 0.5. Data are shown Fig. |l(b)[ they are naturally 
smoother than the true image and the small fluctuations are less visible and corrupted by the noise. 
The empirical mean level of the image is correctly observed (the null frequency coefficient of 

o 

is ho = 1) so the parameter 70 is considered as a nuisance parameter. Consequently it is integrated 
out under a Dirac (see Sec. |4]). This is equivalent to fix its value to in the algorithm Fig. HI line 4. 
Finally, the method is evaluated on two different situations. 

1. The unsupervised and non-myopic case: the parameters w are known. Consequently, there 
is no Metropolis-Hastings step (Sec. l5.Cl) : lines 9 to 16 are ignored in the algorithm of Fig. [9] 
and w is set to its true value. To obtain sufficient law exploration, the algorithm is run until 
the difference between two successive empirical means is less than 10""^. In this case, 921 
samples are necessary and they are computed in approximately 12 seconds on a processor at 
2.66 GHz with Matlab, 

2. The unsupervised and myopic case: all the parameters are estimated. To obtain sufficient 
law exploration, the algorithm is run until the difference between two successive empirical 
means is less than 5 x 10"^. In this case, 18 715 samples are needed and they are computed 
in approximately 7 minutes. 

Remark 2 — The algorithm has also been run for up to 1 000 000 samples, in both cases, without 
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perceptible qualitative changes. 



6.B. Estimation results 
6.B.I. Images 

The two results for the image are given Figs ■ 1 1 (c)] and 1 1 (d)| for the non-myopic and the myopic 



cases, respectively. 

The effect of deconvolution is notable on the image, as well as on the shown profile. The ob- 
ject is correctly positioned, the orders of magnitude are respected and the mean level is correctly 
reconstructed. The image is restored, more details are visible and the profiles are closer matching 
to the true image than data. More precisely, the pixels 20-25 of the 68-th line in Fig. [T] show the 
restoration of the original dynamic whereas it is not visible in the data. Between pixels 70 and 110, 
fluctuations not visible in data are also correctly restored. 

In order to visualize and study the spectral contents of the images, circular average of empirical 
power spectral density is considered and called "spectrum" hereafter. The subjacent spectral vari- 
able is a radial frequency / such as p = + The spectrum of the true object, data and restored 
object are shown Figs. |3(a)] and |3(b)] in non-myopic and myopic cases, respectively. It is clear that 



the spectrum of the true image is correctly retrieved, in both cases, up to the radial frequency 
/ ~ 0.075. Above this frequency, noise is clearly dominant and information about the image is 
almost lost. In other words, the method produces correct spectral equalization in the properly ob- 
served frequency band. The result is expected from a Wiener-Hunt method but the achievement is 
the joint estimation of hyperparameter and instrument parameters in addition to the correct spectral 
equalization. 

Concerning a comparison between non-myopic and myopic cases, there is no visual differences. 
The spectrum Figs. |3(a)] and |3(b)] in non-myopic and myopic cases respectively are visually indis- 



tinguishable. This is also the case when comparing Figs . 1 1 (c)] and 1 1 (d)] and especially 68-th line. 



From a more precise quantitative evaluation, a slight difference is observed and detailed below. 
In order to quantify performances, a normalized euclidean distance 

e = - a;*||/||a;*|| (46) 

between an image x and the true image x* is considered. It is computed between true image and 
estimate images as well as between true image and data. Results are reported in Tab.[T]and confirm 
that the deconvolution is effective with an error of approximately 6 % in myopic case compared 
to 11 % with data. Both non-myopic and myopic deconvolution reduce error by a factor 1.7 with 
respect to the observed data. 

Regarding a comparison between non-myopic and myopic case, the errors are almost the same, 
with a slightly lower value for the non-myopic case, as expected. This difference is coherent with 
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the intuition: more information are injected in the non-myopic case through the true PSF parame- 
ters values. 

6.B.2. Hyperparameters and instrument parameters 

Concerning the other parameters, their estimates are close to the true values and are reported in 
Tab. [2l The 7^ estimate is very close to the true value with % = 0.49 instead of 0.5 in the two 
cases. The error for the PSF parameters are 0.35%, 2.7% and 1.9% for Wa, and ip, respectively. 
The value of 71 is underestimated in the two cases with approximately 1 .7 instead of 2. All the true 
values fall in the /t ± Sa interval. 

In order to deepen the numerical study, the paper evaluates the capability of the method to 
accurately select the best values for hyperparameters and instrument parameters. To this end, 
we compute the estimation error Eq. for a set of "exhaustive" values of the parameters 
[7e, 7i, Wa-, W/^, if]. The protocol is the following: 1) choose a new value for a parameter (7^ for 
example) and fix the other parameters to the value provided by our algorithm, 2) compute the 
Wiener-Hunt solution (Sec. l5.Al) and 3) compute the error index. 

Results are reported in Fig. SI In each case, smooth variation of error is observed when varying 
hyperparameters and instrument parameters and an unique optimum is visible. By this way, one 
can find the value of the parameters that provide the best Wiener-Hunt solution when the true image 
X* is known. It is reported on Tab. [Hand shows almost imperceptible improvement: optimization 
of the parameters (based on the true image x*) allow negligible improvement (smaller than 0.02 % 
as reported in Tab. [T]). 

So, the main conclusion is that, the unsupervised and myopic proposed approach is a relevant 
tool in order to tune parameters: it works (without the knowledge of the true image), as well as an 
optimal approach (based on the knowledge of the true image). 

6.C. A posteriori law characteristics 

This section describes the a posteriori law using histograms, means and variances of the parame- 
ters. The sample histograms. Figs. [5] and [6l provide an approximation of the marginal posterior law 
for each parameter. Tabs.[i|and[2]report the variance for the image and law parameters respectively 
and thus allow to quantify the uncertainty. 

6.C. 1 . Hypeiparameter characteristics 

The histograms for 7^ and 71, Fig. [51 are concentrated around a mean value in both non-myopic and 
myopic cases. The variance for 7^ is lower than the one for 71 and it can be explained as follows. 

The observed data are directly impacted by noise (present at the system output) whereas they are 
indirectly impacted by the object (present at the system input). The convolution system damages 
the object and not the noise: as a consequence, the parameter 7^ (that drives noise law) is more 
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reliably estimated than 71 (that drives object law). 

A second observation is the smaller variance for 71 in the non-myopic case Fig. |5(c)| than in the 
myopic case Fig. |5(d)[ It is the consequence of the addition of information in the non-myopic case 
w.r.t. the myopic one, through the value of the PSF parameters. In the myopic case, the estimates 
are founded on the knowledge of an interval for the values of the instrument parameters, whereas 
in the non-myopic case, the estimates are founded on the true values for the instrument parameters. 

6.C.2. PSF parameter characteristics 

Fig. [6] gives histograms for the three PSF parameters and their appearances are quite different 
from the one for hyperparameters. The histograms for Wa and wp. Figs. |6(a)] and |6(bj] are not as 
concentrated as the one of Fig. [5] for hyperparameters. Their variances are quite large with regards 
to the interval of the prior law. On the contrary, the histogram for the parameter 93, Fig. |6(c)[ has 



the smallest variance. It is analyzed as a consequence of a larger sensitivity of the data w.r.t. the 
parameter (p than w.r.t. the parameters Wa and wp. In an equivalent manner, the observed data are 
more informative about the parameter (p than about the parameters Wa and Wj3. 

6.C.3. Mark of the myopic ambiguity 

Finally, a correlation between parameters (71, Wa) and (71, Wfj ) is visible on their joint histograms 
Fig.|71 It can be interpreted as a consequence of the ambiguity in the primitive myopic deconvolu- 
tion problem, in the following manner: the parameters 71 and w both participate in the interpreta- 
tion of the spectral content of data, 71 as a scale factor and w as a shape factor. An increase of Wa 
or Wj3 results in a decrease of the cutoff frequency of the observation system. In order to explain 
the spectral content of a given data set, the spectrum of the original image must contain more high 



frequencies, i.e., a smaller 71. This is also observed on the histogram illustrated Fig. 7(a) 



6.D. MCMC algorithm characteristics 

Globally, the chains of Figs. [5] and [6l have a Markov feature (correlated) and explore the parameter 
space. They have a burn-in period followed by a stationary state. This characteristic has always 
been observed regardless the initialization. For fixed experimental conditions, the stationary state 
of multiple runs was always around the same value. Considering different initializations, the only 
visible change is on the length of the burn-in period. 

More precisely, the chain of 7^ is concentrated in a small interval, the bum-in period is very 
short (less than 10 samples) and its evolution seems independent of the other parameters. The 
chain of 71 has a larger exploration, the burn-in period is longer (approximately 200 samples) and 
the histogram is larger. This is in accordance with the analysis of Section [6. C. II 

About the PSF parameters, the behaviour is different for {wa, wp) and The chain of the two 
width parameters has a very good exploration with quasi-instantaneous bum-in period. Conversely, 
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the chain of (p is more concentrated and its burn-in period is approximately 4 000 samples. This is 
also in accordance with previous analysis (Section [6 .C. 21) . 

Acceptation rates in the Metropolis-Hastings algorithm are reported in Tab. [3l they are quite 
small, especially for the rotation parameter. This is due to the structure of the implemented algo- 
rithm: an independant Metropolis-Hastings algorithm with the prior law as a proposition law. The 
main advantage of this choice is its simplicity but as a counterpart, a high rejection rate is observed 
due to a large a priori interval for the angle parameter. A future work will be devoted to the design 
of more accurate proposition law. 

6. E. Robustness of prior image model 

Fig. [8] illustrates the proposed method on a more realistic image with heterogeneous spatial struc- 
tures. The original is the Lena image and the data has been obtained with the same Gaussian PSF 
and also corruption by white Gaussian noise. The Fig. |8(b)| shows that the restored image is closer 
to the true one than the data. Smaller structures are visible and edges are sharper, for example 
around pixel 200. The estimated parameters are % = 1.98 while the true value is 7* = 2. Con- 
cerning the PSF parameters, the results are u?^ = 19.3, u?^ = 7.5 and Cp = 1.15 while the true 
values are respectively ly* = 20, w*p = 7 and = 1.05 as in the previous section. Here again, the 
estimated PSF parameters are close to the true values giving a first assessment of the capability of 
the method in a more realistic context. 

7. Conclusion and perspectives 

This paper presents a new global and coherent method for myopic and unsupervised deconvolution 
of relatively smooth images. It is built within a Bayesian framework and a proper extended a 
posteriori law for the PSF parameters, the hyperparameters and the image. The estimate, defined 
as the posterior mean, is computed by means of an MCMC algorithm in less than a few minutes. 

Numerical assessment testifies that the parameters of the PSF and the parameters of the prior 
laws are precisely estimated. In addition, results also demonstrate that the myopic and unsupervised 
deconvolved image is closer to the true image than the data and show true restored high-frequencies 
as well as spatial details. 

The paper focuses on linear invariant model often encountered in astronomy, medical imaging, 
nondestructive testing and especially in optical problems. Non-invariant linear models can also 
be considered in order to address other applications such as spectrometry [l4| or fluorescence mi- 
croscopy [[T3I . The loss of invariance property precludes entirely Fourier-based computations but 
the methodology remains valid and practicable. In particular, it is possible to draw samples of the 
image by means of an optimization algorithm [[381 . 

Gaussian law, related to L2 penalization, is known for possible excessive sharp edges penaliza- 
tion in the restored object. The use of convex L2 — Li penalization [[39] - l4TI or non convex L2 — Lq 
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penalization H2ll can overcome this limitation. In these cases a difficulty occurs in the development 
of myopic and unsupervised deconvolution: the partition function of the prior law for the image 
is in intricate or even unknown dependency w.r.t. the parameters [|Tl|71|43l. However a recent pa- 
per II4TII overcome the difficulty resulting in an efficient unsupervised deconvolution and we plan 
to extend this work for the myopic case. 

Regarding noise, Gaussian likelihood limits robustness to outliers or aberrant data and it is pos- 
sible to appeal to robust law such as Huber penalization in order to bypass the limitation. Never- 
theless, the partition function for the noise law is again difficult or impossible to manage and it is 
possible to resort to the idea proposed in [l4T| to overcome the difficulty. 

Finally, estimation of parameters of correlation matrix (cutoff frequency, attenuation coeffi- 
cients,. . . ) is possible within the same methodological framework. This could be achieved for the 
correlation matrix of the object or the noise. As for the PSF parameters, the approach could rely 
on an extended a posteriori law, including the new parameters and a Metropolis-Hastings sampler. 

8. Acknowledgment 

The authors would like to thank Professor Alain Abergel in IAS laboratory at Universite Paris-Sud 
11, France, for fruitful discussions and constructive suggestions. The authors are also grateful to 
Cornelia Vacar, IMS laboratory, for carefully reading the paper. 

A. Law in Fourier space 

For a Gaussian vector x ~ A/'(/^, S), the law for x = Fx (the fft of x) is also Gaussian whose 
first two moments are the following: 
• The mean is 



A* = E[a;] = FE[x] = F/x. 



(47) 



• The covariance matrix is 



± = E[{x - li){x - liY] = FSFt. 



(48) 



Moreover, if the covariance matrix S is circulant it writes 



S = FSF^ = As. 



(49) 



i.e., the covariance matrix E is diagonal. 



B. The Gamma probability density 



B.A. Definition 



The Gamma pdf for 7 > 0, with given parameter a > and /3 > 0, is written 



6^(7l«,/3) 



1 



.7"-iexp(-7//3). 



(50) 
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Tab.Hlgives three limit cases for (a, (3). The following properties hold: 

• The mean is Eg [7] = a/S 

• The variance is Vg[7] = 

• The maximiser is /3 (a — 1) if and only if a > 1 
B.B. Marginalisation 

First consider a N dimensional zero-mean Gaussian vector with a given precision matrix 7r with 
7 > 0. The pdf reads 

p{x\-i) = (27r)-^/ det[r]^/2 \^_^x'Tx /2] . (51) 

So consider the conjugate pdf for 7 as a Gamma law with parameter (a, (3) (see previous An- 
nex). The joint law for {x, 7) is the product of the pdf given by Eq. (l50l) and Eq. (ISTI) : p{x, 7) = 
p(a; 17)^(7). The marginalization of the joint law is known p4|: 



p{x) = / p{x\-f)p{-f) d7 

_ /3^/Met[r]^/^r(a + iV/2) / /3a;*rx \ """^/^ ^^^^ 
(27r)^/2r(a) 1^ ^ 2 j 

which is a dimensional f-Student law of 2a degrees of freedom with a /3r precision matrix. 
Finally, the conditional law reads: 

P{l\x) = ^^""^"Jf^^'/^^'^' 7"^^/^-^ exp [-7 [x'Tx /2 + 1//?)] . (53) 

Thanks to conjugacy, it is also a Gamma pdf with parameters a , (3 given by a = a + A^/2 and 
= + 2/(a;*ra;). 

C. The Metropolis-Hastings algorithm 

The Metropolis-Hastings algorithm provides samples of a target law f{w) that cannot be directly 
sampled but can be evaluated, at least up to a multiplicative constant. Using the so called "instru- 
ment law" q (^Wp\w^^^^ , samples of the target law are obtained by the following iterations. 

1. Sample a proposition ^ q {w^\w'^^'> 

2. Compute the probability 



3. Take 



^(t+i) ^ J with p probability ^^^^ 

ly*^*) with 1 — p probability . 
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At convergence, the samples follow the target law /{w) [I251I36II . When q (wp\w^^^j = qiwp) the 
algorithm is named independent Metropolis-Hastings. In addition, if the instrument law is uniform, 
the acceptance probability gets simpler in 
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Table 1. Error e (Eq. (1461) ) and averaged standard deviation a of the posterior image 
law. The "Best" error has been obtained with the knowledge of the true image. 





Data 


Non-myopic 


Myopic 


Best 


Error (e) 


11.092 % 


6.241 % 


6.253 % 


6.235 % 


a ofx law 




3.16 


3.25 
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Table 2. Quantitative evaluation: true and estimated values of hyperparameters and 
PSF parameters. 







7e ± 0- 


7i ± 6- 


Wa ± 0" 


± (T 


± (T 




True value 


0.5 


2 


20 


7 


1.05 (7r/3) 


Non-myopic 


Estimate 


0.49 ±0.0056 


1.78 ±0.14 










Error 


2.0% 


11 % 








Myopic 


Estimate 


0.49 ±0.0056 


1.65 ±0.15 


20.07 ±0.53 


7.19 ±0.38 


1.03 ±0.04 




Error 


2.0% 


18 % 


0.35 % 


2.7 % 


1.9% 
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Table 3. Acceptation rate. 



Parameter 








Acceptation rate 


14.50 % 


9.44 % 


2.14% 
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Table 4. Specific laws obtained as limit of the Gamma pdf. 





a 




Jeffreys 





+ 00 


Uniform 


1 


+ 00 


Dirac 
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Fig. 1. The figure [T(a)] represents a 128 x 128 sample of the a priori law for the 
object with 70 = 1 and 71 = 2. Fig. |l(b)| is the data computed with the PSF shown 
in Fig. [2I Figs. |l(c)| and |l(d)| are the estimates with non-myopic and the myopic 
estimate, respectively. Profiles correspond to the 68-th line. 
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Fig. 2. PSF with Wa = 20, = 7 and (f = n/S. The x-axis and y-axis are reduced 
frequency. 




(a) Non-Myopic (b) Myopic 

Fig. 3. Circular average of the empirical power spectral density of the image, the 
convolued image, the data (convolued image corrupted by noise) and the estimates, 
in radial frequency with y-axis in logarithmic scale. The x-axis is the radial fre- 
quency. 
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Fig. 4. Computation of the best parameters in the sense e Eq. (|46l ). The symbol 
' X ' is the minimum and the symbol ' .' is the estimated value by our approach. The 
y-axis of 7^ and 71 are in logarithmic scale. 
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Fig. 5. Histograms and chains for the non-myopic case in Figs. 5(a) 5(c) and the 



myopic case in Figs. |5(b)H5(d)| for 7^ and 71, respectively. The symbol x localizes 
the initial value and the dashed line corresponds to the true value. The x-axis are 
iteration's index for the chains and parameter value for the histograms. 
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Fig. 6. Histogram and chain for the PSF parameters Wa in Fig. [6(a)} Wfi in Fig. |6(b)| 
and Lp in Fig. |6(c)[ The symbol x localizes the initial value and the dashed line 
corresponds to the true value. The x-axis for the histograms and the y-axis of the 
chain are limits of a priori law. 
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Fig. 7. Joint histograms for the couple (71, Wa) and (71, w;^) in Figs. |7 (a)] and |7(b)| 
respectively. The x-axis and y-axis are the parameter value. 
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(a) Data 



(b) Estimated image 



Fig. 8. Observed image Fig. |8(a)| and restored image Fig. |8(b)[ Profiles correspond 
to the 68-th line. The solid line is the true profile. Dashed line correspond to data in 
Fig. |8(a)| and estimated profiles in Fig. |8(b)[ 
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1: Initialisation of x^^\ w^^\ k = 
2; repeat 

% Sample of x 

5: a;^'^ ^ /x + S-^/^ * randn 

% Sample of 7 
6: 7^^'''^ -f- gamrnd(ae, 
7: 7^^ ^ gamrnd(ai,/3i) 

8: 7o^^ ^ gamrnd(ao, /3o) 

% Sample of 

9: lUp ^ rand * {a — b) + a 

10: J ^ 7. (11^ - xf-B- \H,n,, xf) /2 

11: if log(rand) < min{J, 0} then 
12: W^'^^ ^ Wp 

13: \h ^ Aff^^p 

14: else 

15: W^''^ ^ itj(fe-l) 

16: end if 

% Empirical mean 
17: + 1 

18: a;^'VE.^^^Vfc 

19: until 1*^'') — x^''~^^\/\x^''M < criterion 



Fig. 9. Pseudo-code algorithm, gamrnd, rand and randn draw samples of gamma 
variable, uniform variable, and zero-mean unit-variance white complex Gaussian 
vector respectively. 
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